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ABSTRACT 

The origin of spiral patterns in galaxies is still not fully understood. Similar features also develop 
readily in A-body simulations of isolated cool, coUisionless disks, yet even here the mechanism has yet 
to be explained. In this series of papers, I present a detailed study of the origin of spiral activity in 
simulations in the hope that the mechanism that causes the patterns is also responsible for some of these 
features galaxies. 

In this first paper, I use a suite of highly idealized simulations of a linearly stable disk that employ 
increasing numbers of particles. While the amplitudes of initial non-axisymmetric features scale as the 
inverse square- root of the number of particles employed, the final amplitude of the patterns is independent 
of the particle number. I find that the amplitudes of non-axisymmetric disturbances grow in two distinct 
phases: slow growth occurs when the relative overdensity is below ~ 2%, but above this level the 
amplitude rises more rapidly. I show that all features, even of very low amplitude, scatter particles at 
the inner Lindblad resonance, changing the distribution of particles in the disk in such a way as to foster 
continued growth. Stronger scattering by larger amplitude waves provokes a vigorous instability that is 
a true linear mode of the modified disk. 
Subject headings: galaxies: evolution - galaxies: kinematics and dynamics - galaxies: spiral 



1. INTRODUCTION 

Despite considerable effort over many decades, no the- 
ory for the origin of spiral patterns in galaxies has yet 
gained wide accept ance. There is good evidence ( e.g.. Ko- 



t'lat 



rmendy & Norman Tl979l : iKendall et a;.ll201lD 
patterns in galaxies are excited by tidal iiiteract it)ns 
iBvrd fc Howardlll992l:ISalo fc Laurikainenlll993l: 



some 

{e.g., 

Dobbs et 



al. l2010f ). others mav be driven bv bars (e.g.. IS do et al\ 
I2OIOD, and perhaps even d ark matter halo sub itructure 
{e.g., iDubinski et al\ l2008f ) . But when all these agents 
are excluded from A-body simulations, they still manifest 
spiral patterns and it therefore seems reasonable ;o expect 
that self-excited features also develop in galaxies. 

Several reviews (Toomre 1977: A thana s soula 19 g4 : Bcrt in 
& Lin 119961 : iBi'nnev fc Tremaing11J08nSellwoo|Il2010al ) 
give detailed summaries of theoretical efforts to ac :ount for 
self-excited spirals. Most theorists agree that spirals are 
gravitationally-driven density waves in the old stellar disk, 
a point of view that is strongl y supported by the l ight dis- 
trib ution iii the near IR (e.g.. Grosb0l et all 1200' j^Zibetti 
et aL I2OO9I : IKendall et a/.ll201lD and by strea ming m otions 
in high spatial resolution velocity maps (^e.o.. lVisjeiill978l : 
IChemin et Zl l2006l : IShettv et all 120071 ). However, strong 
disagreements arise over other aspects; in particular, the 
lifetimes of the spirals. 

Small-amplitude, non-axisymmetric Jeans instabilities 
of equilibrium models are of two fundamental tyj es. Cav- 
ity modes are standing waves, reminiscent of those in organ 
pipes and guitar strings, in which traveling wavjs reflect 
off an inner boundary, often the galaxy center, and corota- 
tion, where they are amplified. The bar-forming instability 
is the strongest rn ode of many simple models (e.g. , Kalnaja 
119781 : lJalalil2007f ) and is of the cavity type ([ijtoomrc 198$? 
The o ther type of mode , such as the edge mode (,Toomrg 
Il98ll : iPapaloizou fc LinI 119891 ) or groove mode (JLovelace 



fc Hohlfeldl [19781 : iSellwood fc Kahnl[l99Th . is driven from 
corotation. All these are vigorous instabilities that expo- 
ncntiatc or an orbital time scale. 

Bertin & Lin (|1996f ) argue t hat spirals a re caused by a 
more gentk WASER instability ()Markll977f ). a cavity mode 
that reflecis off a "Q-barrier" instead of the center, that 
may becon e quasi-steady at finite amplitude through non- 
li near effec :s. However, their specific examples f Bertin et 
al. 119891 ) O' low-mass, cool disks that are needed to sup- 
port slowl y -grow ing bi-symmetric modes have been shown 
(|Sellwoodl [ ^011h to develop more vigorous multi-arm in- 
stabilities that quickly alter the b asic s tate they invok e. 

Goldreicn & Lynden-BcU (|1965[ ) and lToomrel (|l990t ) ar- 
gue that, i:i galaxies lacking a bar or a recent interaction 
with a coir panion, spirals are no more than shearing bits 
and pieces of swing-amplified noise. Substantial inhomo- 
geneities, s ich as star clusters and giant molecular clouds, 
raise the nc lise level far above that expected from the stars 
alone, and the noise level can be amplified ^ 100-fold in 
cool, responsive disks. However, since the spirals they wish 
to account for themselves cause the disk to heat quickly, it 
is hard to i nderstand how a galaxy could sustain the com- 
bination of responsiveness and noise level needed to create 
spiral featr res of the amplitude we observe. Furthermore, 
this mecha lism has difficulty accounting for even a mod- 
est degree of svmmetrv in the patterns (e.g.. Kendall et 

«'• 120111). 

From th ! earliest davs (iM iller e ^ allll971l: Hocknev fc 
Brownrigg" 19741 : iJames fc S ellwood Il978f ). A-bodv sim- 
ulations of dynamically cool coUisionless disks have man- 
ifested, recurrent short-lived transient spirals. Claims of 
long-lived ' raves in some simulations are not reproducible 
(ISellwoodl E ioTl . ISellwood fc Carlberd (119841 ) showed that 
spiral acti'^ity in purely stellar disks is self-limiting, be- 
cause potei itial fluctuations caused by the rapidly evolving 



spirals themselves scatter stars away frora circular orbits; 
the increased random motion in the disk reduces its ability 
to support furth er collective oscillations. Sellwood & Carl- 
berg y,984) also found that spiral activity can continue 
"indefinitely" in models that mimic the effects of gas dis- 
sipation and star formation, offering an attractive expla- 
nation for the observed coincidence between spiral activity 
and gas content in galaxies. The behavior they reported 
continues to occur in modern simulations of gala xy forma- 
ti on (e.q..lAbadi et aLll2003l:lRoskar et a/.ll200a Agertz et 
aL[2011f). 

Local theory ([Julian fc Toomrg Il966f ) predicts that a 
perturbing mass moving on a circular orbit in a self-gravitating 
disk surrounds itself with a wake. In an A^-body simula- 
tion, each particle is itself a perturber and has its own 
peculiar motion that is typically as large as that of the 
surrounding particles whose motion it affects. Note that a 
system of dressed particles, with each orbiting at its own 
angular rate, will give rise to shearing density fluctuations 
that will be strongest when the wakes are aligned, giving 
the appearance of swing amplification. In fact, the picture 
of randomly placed dressed particles in a shearing disk is 
equivalent to that of swing-amplified noise, with the in- 
put noise level raised by the m ass cloud surrounding each 
source particle (iToomrd 1990t) . 

Toomre & Kalnajs ( 199lf ) studied swing-amplified par- 
ticle noise using simulations in shearing (x, y)-coordinates 
in a periodic box. They had to counter the continuous 
rise in the level of random motion in order to be able to 
study the long-term equilibrium level of spiral activity, and 
they therefore added a "dashpot" type damping term to 
the radial acceleration, f^ — —Cvx- The level of random 
motion should have settled to a roughly constant value if 
the damping constant C were proportional to the inverse 
of the particle number. They found, however, that the 
observed rms velocity of the particles exceeded their pre- 
dicted level, because the density fluctuations were about 
twice as large as they expected - a discrepancy that they 
attributed to additional correlations between the particles 
that developed over a long period. 

While the origin of spiral features in global simulations 
has yet to be explained, their amplitude in larger A-body 
simulations see ms to be too great to be just s wi ng- amplified 
artic le noise (" Sellwood fc CarlberS Il984l ). iFuiii et all 
20111 ) found that spirals in simulations with larger num- 
bers of particles heated the disk more slowly, allowing the 
spiral activity to continue for longer. These authors incor- 
rectly attribut ed the rapid heating rate in t he smaller A 
simulations by ISellwood fc Carlberel (|1984[ ) to two-body 
relaxation, whereas the higher level of shot noise simply 
provokes stronger spirals Q 

In this paper, I attempt to uncover the reason for the 
suprisingiy vigorous spiral activity in A^-body simulations 
of cool, coUisionless disks. The models employed here are 

^ If coUisional relaxation were rapid, it would have been impossible 
for the same codet^ia^^uggortedthc global modes (Sellwood & : 
Athanassoula 19861: lEarn fc Sellwoodl ll996: Sellwood k. Evans 2001D 
predicted by linear perturbation theory for a coUisionless stellar fluid, 
or to match the results fro m direct integ ration of the coUisionless 
Boltzmann equation (Inaga ki et aZ.iri986D . The formula for the re- 
laxation rate in a two-dimensional disk of particles does not contain 
the Coulomb logarithm (Rybicki 1972) and the rate is therefore sup- 
pressed more effectively by gravity softening. 



highly idealized in order to simplify the dynamics to the 
point that it becomes understandable. Improved realism 
is deferred to later pa pers. 

I have long hoped (|Sellwoodl[T99ll . |2000H that, in real- 
istic galaxy models, the decay of one spiral feature might 
change the background disk in such a manner as to cre- 
ate coiiditioiis_Jo£^^^ie^^nstabilit^j_^nu^^ & 
Lin (1989) had observed in simulations of a cool, low-mass 
mass disk having a near-Keplerian rotation curve. Here I 
finally present the evidence that something similar to this 
idea really does occur in more galaxy-like models, but the 
behavior differs in detail from that I envisioned. 

2. A LINEARLY STABLE DISK 

The disk in which the circular speed is independent 
of radius , now known as the "Mestel disk" ( Binnev & 
Trema ine 20(381') . is predicted by linear stability a nalyses 
(|Zangj 119761 : iToomTi Il98ll : lEvans fc ReadI Il998[ ) to be 
globally stable to most non-axisymmetric disturbances. I 
therefore use this model in the simulations described in 
this paper. 

The potential in the disk plane is 



$(i?) = K,^ln(-| 



(1) 



where Vo is the circular speed, and Ri is some reference 
radius. The surface density of the full-mass disk is 



S(i?) = 



1^2 



(2) 



2ttGR 

A star at radius R moving in the plane with velocity com- 
ponents vr and v^ in the radial and azimuthal directions 
respectively, has specific energy E = $(i?) + {vj^ + f5)/2 
and specific angular momentum Lz = Rv,j,. A distribution 
function (DF) for this disk that yields a Gaussian distri- 
bution of radial velocitie s of width a^ is (|Tbomrdll977l : 
iBinnev fc Tremaind[2008h 



fiE,L) 



otherwise, 



(3) 



where q = Vq /(^r ^ 1 and the normalization constant 
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Vi/i2nG) 



2«/V^(9/2-0.5)H+'' 



(4) 



In order to avoid includin g stars havin g very high orbital 
frequencies near the center. IZangl (|l976l ) introduced a cen- 
tral cut out in the active surface density by multiplying 
the distribution function / by the taper function 

^'^^'^ ^ Jr^VoTT^' ^^^ 

where the index v determines the sharpness of the taper. 
This taper function removes low angular momentum stars, 
reducing the active disk mass at the center, which must 
be replaced by unresponsive matter in order that the total 
potential is unchanged. The characteristic radius Ri of 
this cut out introduces a scale length. Henceforth, I adopt 
units such tha t Vq = Ri — G — 1. 

Zang (1970) found that this almost full-mass disk was 
stable to all bisymmetric modes, provided the index v < 2. 
This important result was confirmed by lEvans fc ReadI 
(|1998r ). who also found that other disks with power-law 
rotation curves were remarkably stable to bisymmetric 



modes. We now understand that this disk model is linearly 
stable to bar-forming modes through Toomre's (1981) mech- 
anism, because the frequency singularity at the disk cen- 
ter ensures that disturbances of any pattern speed flp will 
have an inner Lindblad resonance (hereafter ILR) as long 
as 771 > 1. 

However, Zang also found that the full-mass Mestel disk 
was globally unstabl e to lop-sided (m = 1) modes, as are 
the power-law disks (lEvans fc Readl[l998| ) . This instability 
persists no matter how hot the disk or gentle the central 
cut out. The mechanism for the lop-sided mode is also a 
cavity mode, since feedback through the center is allowed 
only for m = 1 waves because the ILR condition in power- 
law disks cannot be satisfied when Xlp > 0. 

Swing-amplification in a disk ha ving a flat rotation curve 
is most vigorous for 1 < X <2.5 (jJulian fc Toomrg|1966l : 
iTbo mrc 1981), where X — Rkcrit/'m, with the smallest 
wavenumber for axisymmetric instabilitie s fccrit — K^/(27rGE) 
(|Toomrdll96"l : lBinnev fc Tremainell2008f ). Since X = 2/m 
in a full-mass Mestel disk, the swing-amplifier is at its 
most vigorous for m = 1 disturbances. The parameter 
X = 2/{^m) is increased by reducing the active surface 
density of the disk to ^E without changing the circular 
speed; here < ^ < 1. Since the swing-amplifier is ail-but 
dead when X > 3 in a flat rotation curve, the lop-sided 
instability should disappear for ^ < 2/3. 

I'oomre Il981 ) reports that a half-mass (^ — 0.5) disk, 
with v — 4: for the central cut out, has no small-amplitude, 
global instabilities for any m. This model is the only 
known globally stable disk with differential rotation and 
significant disk massQ 

As th e disk still extends indefinitely to large R, iToomrd 
(|cl989l ) added an outer taper 



To{L,) = 



1 



RoVo 



-\ -1 



(6) 



with index fi to control the abruptness of the taper about 
the angular momentum characteristic of a circular orbit at 
some large radius Ro- Again, the potential is unchanged 
if the mass removed is replaced by rigid matter. He found 
At < 6 to be adequate to avoid provokin g outer edge modes 
(|Toomrdll98ll : iPapaloizou fc Linlll989f ) . 

As this taper function still does not limit the radial ex- 
tent of the disk, but merely reduces the active density of 
the outer disk, I have to add the additional constraint that 
f{E,Lz) = when E > $(i?max), so that the disk con- 
tains no particles with orbits that extend beyond -Rmax- 
In order that this additional limit does not also provoke 
an edge mode, I choose i?max ^ 1.7i?o, so that it removes 
mass only where the disk density is already very low. 

In the simulations reported in this paper, I adopt ^ = 
0.5, q = 11 AA so that the half- mass disk has nominal Q — 
1.5, the inner cut out is characterized hy ly — A at Ri ~ 1, 
and the outer taper has /i = 5 at i?o = 11. 5 and the 
limiting radius is -Rmax = 20. 

Figure [T] shows the de nsity of particles in the spac e of the 
two actions Jr and J^ (jBinney fc Tremaindl2008l ). In an 
axisymmetric two-dimensional disk model, the azimuthal 



^ Two of the composite fl models reported bv lKalnaisI II1972I') appear 
to be the only other known stable disks, but are of less interest 
because they rotate uniformly. 




Fig. 1 . — Initial density of particles in the space of the two actions. 
The contour spacings are linear from 5% to 95% of the function 
maximum. The density for low values of J^ is suppressed by the 
inner angular momentum taper and, as is generally true, the density 
is greatest for near circular orbits ( J^ = 0) and decreases for more 
eccentric orbits. 



action J^ — L^, while the radial action, 

Jr = - vr dR, (7) 

expresses the amount of radial motion the particle pos- 
sesses. Here vr = [2E - 2$(i?) - Ll/R^] ^^^ (positive 
root only in this definition) and the limits of integration 
are the radii of peri- and apo-center of the orbit where the 
argument of the square root is equal to zero. Thus Jr also 
has the dimensions of angular momentum and Jr = for 
a circular orbit. The actions are an alternative set of inte- 
grals to the classical integrals E and L^, but the relation 
between the two sets of integrals is cumbersome. 

3. SIMULATIONS OF THE "STABLE" DISK 

The model just described, with the outer tap er but not 
the energy cut-off, was predicted by iToomrd (|cl989) to 
have no small-amplitude instabilities. Here I report sim- 
ulations with this model both to test the prediction and 
to study how replacing the infinitely finely divided "stellar 
fluid" by a finite number of particles affects the behavior. 

I use the two-dimensional pola r-grid code ( Scl lwood 1981) 
with multiple time step zones ()Sellwoodlll985l ). that has 
previously been shown to reproduce the predicted normal 
modes of unstable disks (ISellwood fc Athanassoulal 1 19861 : 
lEarn fc Sellwoodl[T996l : ISellwood fc Evansll200in . "Cravita- 
tional forces between particles obey a Plummer softening 
law, with the softening parameter e = Ri/8, so that forces 
can be thought of arising from a disk of thickness e. 

The grid has 128 spokes and 106 radial rings, and I re- 
strict disturbance forces to the m = 2 sectoral harmonic 
only; the axisymmetric radial force is —Vq/R at all radii. 
I employ a basic time step of St = Ri /{AOVq), which is 
increased by successive factors of 2 for particles with radii 
R > Ri, R > 2Ri, R > 4i?j, k R > 8Ri. Timesteps are 
also sub-divided when pa rticles move near the ce nter in 
a series of "guard zones" IShen fc SellwoodI (|2004f ). with- 
out redetermining the gravitational field, which is an ad- 
equate approximation since the dominant force in this re- 
gion arises from the rigid component. Direct tests showed 
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Fig. 2. — Upper panel: Time evolution of the peak overdensity 
in a series of simulations of the half-mass Mestel d isk with d ifferent 
numbers of particles. The model is predicted by IToomrd l lcl989l) 
to be globally stable. The ordinate reports the maximum value of 
(5S/S on grid rings over the range 1.2 < R < 12 where the surface 
density is little affected by the tapers. Linear theory predicts the 
amplitude should remain proportional to N~^'^. 
Lower panel shows the radius of the measurement reported in the 
upper panel for run 500M only. The apparent lower bound is because 
values for R < 1.2 were excluded to eliminate shot noise in the 
tapered part of the disk. 



and, since I have chosen units such that Ri — Vq — 1, the 
orbit period at radius R is 2TrR. 

As more particles are employed, the initial value of d^ax 
decreases as A^~^/^, as it must for randomly placed par- 
ticles. Its value rises by a factor of a few in the first few 
dynamical times as the initial shot noise is swing-amplified. 
In the smallest N experiment, the initially amplified noise 
is already close to the saturation level, but as N rises, the 
amplitude takes increasingly long to reach its maximum. 
However, the final amplitude is independent of N for all 
particle numbers up to A^ = 5 x 10^. In the larger N exper- 
iments, the amplitude increase is characterized by roughly 
exponential growth at two distinct rates: an initial period 
of slow growth, followed by a steeper rise once Jmax ^ 0.02. 
These two separate rates of growth are approximately in- 
dependent of N over the amplitude ranges where they are 
observed. 

It might seem that exponential growth indicates an un- 
stable normal mode, or perhaps a few such instabilities, 
but such an interpretation is unattractive for a number of 
reasons. Simulations of linearlv unstable models (e.g.. Sell- 
wood & Athanassoula il986 : Earn &: Sellwood 19961 ) usu- 
ally reveal the mostly vigorous instability emerging from 
the noise at an early stage and dominating the subsequent 
growth until it saturates. Modulated growth could occur 
in a system that supports a small number of overstabilities 
having similar growth rates but different pattern speeds, 
since the changing relative phases of the modes over time 
causes beat-like behavior. In all such cases, however, the 
overstabilities should maintain phase coherence until they 
saturate, but power spectra of these simulations (shown 
below) reveal multiple features none of which retains phase 
coherence throughout the period of growth. Furthermore, 
were the later, more rapid rise due to the emergence of 
more rapidly growing normal modes of the original disk, it 
is unusual that they should take so long to rise above the 
amplitude of the more slowly growing "modes" that would 
also need to be invoked to account for the initial slow rise, 
and it is most unlikely that the change of slope would oc- 
cur at almost the same amplitude, as suggested for the 
three larger A^ experiments (Figure [5]). Finally, Toomre's 
(1981) linear stability analysis that predicted the model to 
be globally stable also argues against this interpretation of 
the results. 



that the results are insensitive to moderate changes of grid 
size or time step. 

This paper reports results from several simulations that 
differ in the number of particles, and I label each run by 
the number of particles. Much of the analysis focuses on 
run 50M, which therefore has 5 x 10^ particles, and its 
variants runs 50Ma, 50Mb, etc. 

4. RESULTS 

Figure [D shows the evolution of the relative overdensity 
of bi-symmetric disturbances in a series of simulations with 
increasing numbers of particles. The upper panel shows 
the quantity ^max, which is the largest value of the ratio 
^E/E over the radial range 1.2 < i? < 12, i.e., the part 
of the disk where the surface density has not been signif- 
icantly reduced by the tapers. The unit of time is Ri/Va 



4.1. Slow growth 

The lower panel of Figure [2] shows that the maximum 
disturbance density in run 500M is generally in the inner 
disk for t < 2500, i.e., until the saturation amplitude is 
reached. The behavior over the entire early period < 
t < 2000 seems to follow a quasi-repetitive pattern: a 
density maximum appears in the range 3 < i? < 6 that 
then propagates inwards. This pattern results from swing- 
amplified noise creating a trailing spiral disturbance that 
travels inwards at th e group velocity , as shown in the "dust 
to ash es" Figure of iToomrg ()l98lD . ISalo fc LaurikainenI 
(J2000l , their Figure 7) report very similar behavior in their 
three-dimensional simulations of a low-mass exponential 
disk embedded in a rigid halo. 

The gradual rise of spiral activity in the global simu- 
lations reported here can also be viewed as the build-up 
of mass clouds surrounding each particle. Figure |3] shows 
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Fig. 3. — Gradual development of polarization in the mass distribution of run 50M. The contours show the positive part only of the 
disturbance density, with the bisymmetry enforced by the code, and contour levels are the same in each panel. The inner boundary is at 
R = 1.5, the outer at R = 15, and the "source" particles are at _R = 7. The peak relative ovcrdcnsity at each time is reported in the lower 
right corner. 
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Fig. 4. — Power spectra of m = 2 disturbances in run 50M. Successive contour levels, which start from the same value in each plot, differ 
by factors of two. Each panel is taken from data over the period indicated, which overlap so that only half the plots are from independent 
data. The solid curves show 2Q{R) and the dashed curves 2Q{R) ± k{R). Data outside the radial ranges marked by vertical dotted lines, 
where the surface density is low, are excluded because they are too noisy. 



that each particle can be regarded as being "dressed" by a 
spiral wake, as in the local context. To construct this Fig- 
ure at each instant, I stacked copies of the grid-estimated 
density distribution from run 50M by scaling the entire 
grid radially and rotating it such that the density maxi- 
mum on each "source" ring in turn lay at i? = 7 and zero 
azimuth. At the initial moment, the non-axisymmetric 
part of the combined density shows no features other than 
the azimuthal extension of the source caused by forcing 
cos 29 angular dependence. As time progresses, each source 
particle becomes dressed by a wake of gradually increasing 
mass and spatial scale; the maximum over-density shows a 
general rise, with fluctuations that closely follow the time 
variations of 5max in the same simulation, shown by the 
cyan line in Figure [5] 

However, two related aspects of the behavior are espe- 
cially noteworthy. First, growth in the large N simulations 
continues for a niuch longer period than expected. Julian 
& Toomre (|l966f ) found, from a local linear analvsis. that 
the wake takes ~ 5 epicycle periods to become fully de- 
veloped. Swing-amplification of particle shot noise causes 
the immediate rise in the first few time units, as already 
noted, which is a large part of the development of spiral 
wakes. The epicyclic period at radius R in the Mestel disk 
is 2^/^7r i?/Vb, or ~ 44 time units at i? = 10, about halfway 
out in the disk. Thus the period of continued gradual 
growth greatly exceeds five epicycle periods at a typical 
radius, which is the timescale expected from linear theory. 
Second, the simulations do not reveal a limiting amplitude; 
instead more rapid growth takes over in every case, even 
when 5 x 10* particles are employed. Section l476l presents 
further analysis of the slow-rise phase that accounts for 
these differences. 

4.2. Rapid rise 

The acceleration in the rate of growth once (Jmax ^ 0.02 
is clearly inconsistent with the linear theory prediction. 
The enhanced growth rate is again approximately inde- 
pendent of N, as is the final amplitude at a relative over- 
density between 20% and 30%. 

Figure |4] shows a number of power spectra, i.e., contours 
of power as functions of radius and frequency, from run 
50M. A horizontal ridge in these figures indicates a coher- 
ent density disturbance that extends over a range of radii 
and has an angular frequency mQp, with m — 2 and Qp be- 
ing the rotation rate or pattern speed. There is very little 
power at frequencies higher than those shown. The solid 
line indicates the frequency of circular motion, mfl and 
the dashed lines mil ± k. Each panel gives the spectrum 
over an interval of 200 time units, with the start times of 
successive panels shifted forward by 100 time units. 

As shown in the lower panel of Figure [2] for run 500M, 
density fluctuations at early stages in run 50M also arise at 
radii 4 < i? < 7 from swing-amplified shot noise and sub- 
seque ntly propagate inward at the group velocity (jToomrd 
119691 ). The amplitude of each disturbance rises as it ap- 
proaches the ILR because wave action is conserved as the 
group-velocity decreas es until wave-parti cle interactions 
absorb it near the ILR ()Lvnden-Bell fc Kaln ais 1972; Mark 
IT974h . 

The behavior changes around t ^ 1300. The strongest 
feature up to this time has a frequency ~ 0.5 {e.g., third 
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Fig. 5. — Time evolution of the peak overdensity in five simula- 
tions with N = 50M. The cyan curve is for run 50M, reproduced 
from Figure [2] while the green, red, and magenta curves show the 
behavior of runs 50Ma, 50Mb, and 50Mc created by randomizing the 
azimuths of the particles at t = 1000, 1200 and 1400 respectively. 
For the blue curve, run 50Md, the particles preserved the same E 
and Lz, but the radial phase was also reset. Both runs restarted 
from t = 1400 manifest a single dominant instability. 



panel of Figure 2]), but thereafter a new much stronger 
disturbance appears with a frequency ~ 0.55. Analysis 
with the mo de fitting apparatus described in Sellwood & 
Athanassoula (|l986f ) using data from the period 1100 < 
t < 1500 reveals that this wave has a frequency ui = 
mi^p + i'j ~ 0.55 -1-0. 13i, although the co-existence of other 
waves over the same period makes this frequency estimate 
somewhat uncertain. The significant growth rate suggests 
it is a true instability, yet it is scarcely conceivable that 
this mode could have been growing at this rate from t = 0, 
since that would require its initial amplitude to have been 
smaller by > 50 e-folds before it first became detectable. 

4.3. Randomized restarts 

Figure [5] reports the results of four further simulations 
that were restarted from rearrangements of the particle 
distribution at selected times from run 50M. The original 
run is reproduced as the cyan curve. 

To construct three of these simulations, I simply added 
to the azimuthal phase, a random angle selected uniformly 
from to 27r, while leaving the radius and velocity com- 
ponents in cylindrical polar coordinates unchanged. This 
procedure erases any non-axisymmetric feature above the 
shot noise level, and the starting amplitude is similar to 
that of the cyan line at t = 0. 

The green, red and magenta curves respectively show 
the evolution of runs 50Ma, 50Mb, and 50Mc, which were 
re-started from t = 1000, t = 1200, and t = 1400. In 
all cases the peak non-axisymmetric density rises more 
quickly than that of the original run over the period 200 < 
t < 1000. Furthermore, the later the restart time, the 
more rapidly the amplitude rises and the sooner the am- 
plitude saturates! The different behavior of these random- 
ized runs can only be a consequence of the changes to the 
particle distribution that occurred in the original run prior 




Fig. 6. — Best fit mode extracted from the period of exponential 
growth in rur^50Mcjasing_^he_^rocedure_^escribe^^^Sellwood & 
Athanassoula l|1986f) . The circles mark the radii of the principal 



resonances. 



to azimuthal randomization. 

The most dramatically different behavior is that of run 
50Mc (magenta), which manifests a single vigorous insta- 
bility of eigenfrequency to = 0.6-|-0.16i - close to that of 
the dominant disturbance around t = 1400 identified in 
the earlier analysis of run 50M (Section 14.21) . The form of 
the unstable mode is illustrated in Fig|6l Clearly, the par- 
ticle distribution at t = 1400 has a feature that provokes 
this instability, and which must have been created by the 
earlier evolution. 

The behavior of the other two cases (green and red 
curves of Figure [5]) differs not quite so markedly from the 
original, but again indicates that the particle distribution 
at both t — 1000 and t = 1200 contains features that 
lead to an accelerated rise of non-axisymmetric features. I 
return to these cases in Section WM 

Note that randomizing the azimuthal phases at i = 
1400 only made the identifiable instability of the origi- 
nal model stand out more clearly. Thus, the destabilizing 
agent cannot be a previously-created non-axisymmetric 
feature, e.g., a weak bar, in the particle distribution, and 
must therefore be a feature in the distribution either of 
the radial phases, or of the integrals {E and Lz or Jn and 

In order to test whether changes to the distribution of 
radial phases are important, I determined E and L^ for 
each particle and then selected at random a new pair of 
radial and azimuthal phases from an orbit having these 
integrals in the analytic potential of the Mestel disk. These 
new phases determine new positions and velocities for each 
particle, although it has the same integrals as before. The 
evolution in this case, run 50Md shown by the blue line in 
Figureini makes it clear that the same instability is present 
as that in run 50Mc (magenta line). Thus the important 
change by t = 1400 of the original run is to the distribution 
of the integrals. 




Fig. 7. — Action-space density of particles in run 50M at t = 1200 
(upper) and t = 1400 (lower). 



4.4. Changes to the distribution function 

Figure [7] shows the distribution of particles in action 
space of run 50M at the times t = 1200 (upper) and t = 
1400 (lower). Although the contours are noisier, caused 
by a degradation of the original smooth arrangement, their 
overall shape in the upper panel is little changed from that 
at t = (Figure [T]). However, more substantial changes 
are visible in the lower panel. The changes are shown 
more clearly in the upper panel of Figure |51 which con- 
tours differences between the distributions at i = 1400 
and at i = 0, with increases shown by the blue contours 
and decreases by red. 

The various lines in Figure [5] show resonance loci and 
scattering trajectories for a disturbance of angular fre- 
quency mflp = 0.5, which is the frequency of the dom- 
inant feature of the power spectrum shown in the third 
panel of Figure ID The resonance lines (dashed) are com- 
puted for flp = flfj, + Iflfi/m, where I = for corotation 
and I ~ ^1 for the inner and outer Lindblad resonances, 
and the frequencies fl^ and flu are computed for orbits 
of arbitrary eccentricity. The scattering trajectories (solid 
lines) are computed assuming AE = flpALz, as required 
by the conservation of the Jacobi constant in a rotatin g 
non-axisymmetric potential iBinnev fc Tremaind ()2008[ ). 
with AE converted to AJu and assuming J^, = ini- 
tially for a particle in resonance. Note that any scattering 
at corotation would occur with AJu — 0. 
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Fig. 8. — Upper panel contours the differences between the action- 
space density of particles in run 50M at t = 1400 and t = 0, with 
positive (negative) differences being indicated by blue (red) con- 
tours. The dashed lines show the loci of Lindblad resonances (cyan) 
and of corotation (green) for f2p = 0.25, while the solid lines indicate 
the scattering trajectory from J^ = for the Lindblad resonances, 
lower panel shows the mean displacements of particles in the same 
space and over the same period. The arrow size is proportional to 
the mean value of the displacement in each element. 



It is clear that the principal changes to the distribution 
of particles in action space were caused by scattering at 
the ILR of the strongest non-axisymmetric feature to have 
developed before t = 1400. This conclusion is consider- 
ably strengthened by the evidence in the lower panel of 
Figure [U which shows the mean displacements of particles 
in action space between t — and t — 1400. This Fig- 
ure not only confirms the convergence of vectors towards 
the density maximum created by scattering at the ILR, but 
also reveals that substantial changes have taken place both 
at corotation and at the OLR that produce no net change 
in density. Note that the vectors along the Lindblad res- 
onance lines are parallel to the scattering directions for 
those resonances, computed for the same pattern speed, 
while they are near horizontal at corotation. Thus parti- 
cles at all three resonances are scattered with no change 
to Jacobi's invariant for the measured pattern speed. The 
displacements at the OLR cause no apparent net change 
to the density (upper panel), while those at the ILR do. 
The different behavior at the two Lindblad resonances is 
partly caused by the geometric consequence of the "fo- 
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Fig. 9. — Second moment of the radial velocity distribution of 
particles in run 50M at the three indicated moments. No visible 
changes have occurred by t = 700, while a significant increase can 
be seen near the center at t = 1400. 



cusing" of waves that propagate inwards, while outgoing 
waves spread out over an increasing disk area. It is also 
significant that the scattering direction at the ILR is very 
close to the resonance locus, so that particles stay on res- 
onance as they lose angular momentum, whereas a small 
gain or loss of angular momentum at the OLR moves the 
particle off resonance. 

There are a few additional features visible in the lower 
panel of Figure [5] that appear to be due to another some- 
what lower frequency wave that is also present at i = 1400 
(see the fourth panel of Figure U]) . These features are much 
weaker in an equivalent plot (not included here) of the 
lower panel of Figure |8] constructed for t = 1300. 

Scattering at the ILR naturally causes a very localized 
increase in the random motions of the particles. Figure IHl 
shows the radial variation of the second moment of the ra- 
dial velocities estimated from the particles. By construc- 
tion from the analytic DF the radial velocity distribution 
is near Gaussian ed t — 0, but the resonant scattering in 
action space does not preserve this property. Note that 
this very localized heating at a resonance was caused by a 
disturbance that had a relative overdensity of ^ 2%. 

4.5. Nature of the new instability 

What causes the new instability that appears in runs 
50Mc and 50Md, and which is also present, but less dis- 
tinct, in the continued run 50M? Since randomizing the az- 
imuthal coordinates makes the instability stand out more 
clearly, it cannot be caused by trapping of particles into 
some non-axisymmetric feature, neither can it be caused 
by the kind of non-linear wave coupling discussed bv Mas- 
set & Tagger (19971) or |Fuclis_eLa/J (|2005D. Its rapid expo- 
nential growth characterizes it as a vigorous, global, linear 
mode of the modified disk. ^^_^^^ 

It has long been my expectation {e.g., ISellwoodlll991l . 
l2000t ) that the rec urrent cycle w ould prove to be one of 
groove-type modes (jSellwood fc K ahn 1991) caused by de- 
populating the DF near the OLR of one wave to provoke 
a new instability wit h corotation nea r the p revious OLR. 
Indeed this is what ISellwood fc LinI ([19891 ) observed in 
their simulations of a low mass disk with a near Keplerian 
rotation curve. Yet I have been unable to find evidence 



to support this expectation in more massive disks with 
flatter rotation curves, and the new evidence presented in 
Figure[5]seems quite emphatically to rule it out, at least for 
this case. While changes to the DF are clearly responsible, 
it is not scattering at the OLR, and some other mechanism 
is at work. 

The only significant change to have occurred to cause 
runs 50Mc and 50Md to behave so differently from the 
initial behavior of run 50M is the localized heating near 
R = 1.17, which is the radius of the ILR for near circular 
orbits for the scattering wave. There are no net changes 
elsewhere, even though the lower panel of Figure [S] indi- 
cates that some particles in other parts of the disk have 
changed places. 

Since the new instability has a slightly higher pattern 
speed than that of the preceding scattering wave, the radii 
of the principal resonances are all correspondingly smaller, 
and therefore the ILR of the new instability would lie 
within the heated region of the disk. The localized de- 
ficiency of nearly circular orbits created by the previous 
changes to the DF will have strongly inhibited a support- 
ing response in this region of the disk to forcing by an 
inward propagating density wave. 

Thus it seems to me that the most plausible unstable 
mode created by these "initial" conditions is a cavity-type 
mode that is a standing wave between corotation and a 
hard inner reflection off the region of the disk where the DF 
has been changed. Since the inner disk is "hotter" , there is 
a superficial resemblance between this suggestion and the 
WASER mechanism p r opose d by Mark (1977) and advo- 
cated bv lBertin fc LinI (|l996f ). However, the WASER mode 
has a low growth rate because the out- and in-going waves 
at corotation are both trailing, resulting in very mild am- 
plification as the disturbance connects only the long- and 
short-wave branches of the dispersion relation. Further- 
more, their mode requires a low surface mass density and 
Q ^ 1.0 near corotation, whereas this model has a massive 
disk with Q — 1.5. Finally, the change in disk properties 
that causes the inner reflection is quite different from the 
simple increase in Q postulated for the WASER mecha- 
nism. Much more vigorous modes, such as that observed 
here, se em likely to re quire a full swing from leading to 
trailing (JToomrelll98l[ ). Thus I suspect that the in-going 
trailing wave reflects off the region where the DF has been 
changed as an outgoing leading wave. I have no evidence 
that a hard reflection can occur in this case, so this pro- 
posed mechanism is entire speculation. 

4.6. Earlier restarts 

Simulations 50Ma and 50Mb were restarted at i = 1 000 
and t — 1 200 respectively. As shown in Figure [SJ density 
fluctuations in these runs grew at rates that were interme- 
diate between that of the original run 50M and the clear 
instability exhibited by those (50Mc and 50Md) that were 
restarted at t = 1400. Moreover, these two runs support 
multiple disturbances having several different angular fre- 
quencies, none of which maintains coherence through to 
the saturation amplitude. The randomized particle dis- 
tributions at these intermediate times must also differ in 
important ways from the smooth particle distribution of 
run 50M at t = 0. 

Figure [TU] shows that the changes to the DF hy t = 1000 
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Fig. 10. — Upper panel contours the differences between the 
action-space density of particles in run 50M at i = 1000 and t = 0, 
with positive (negative) differences being indicated by blue (red) 
contours. The kernel width is five times larger than that used in 
Figure [8] and the contour levels are one fifth. 

Lower panel shows the radial velocity dispersion of particles in run 
50M at the three indicated moments. The changes from t = 
are non-zero only near the center; they are tiny by t = 1000 and 
moderately larger at t = 1200. 



are tiny, but nevertheless significant. I had to reduce the 
contour levels in order to reveal the tiny changes to the 
density of particles in action space (top panel) , which ne- 
cessitated a proportionate increase in the kernel size so 
that the shot noise from the number of particles contribut- 
ing to the lowest contour is unchanged. It is clear that 
there has been some significant scattering of particles away 
from the Jr — axis, especially at small radii. These, and 
the significantly larger changes that occur by i = 1200 
(not shown) caused some heating of the innermost part of 
the disk, as illustrated in the lower panel. This very mild 
heating of the inner disk is caused by the absorption at the 
ILRs of incoming disturbances created by swing- amplified 
shot noise at larger radii. Although the wave amplitudes 
were tiny, they caused a lasting and significant change to 
the DF. 

The randomized restarts at i = 1000 (run 50Ma) and 
t ~ 1200 (run 50Mb) confirm that these tiny changes, 
and not non-linear coupling say, were responsible for the 
enhanced growth of non-axisymmetric structure. It seems 
likely that the small changes that have taken place are suf- 
ficient to cause some partial reflection of waves incident on 
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the inner disk that boosts the density fluctuations, but is 
apparently insufficient to provoke an indefinitely growing, 
more vigorous mode at these times. 

Thus the slow rise in the density fluctuations described 
in Section HTTJ occurred because the physical system gradu- 
ally develops a marginally modified inner disk. Mild scat- 
tering of particles in the inner disk as the wave action of 
the noise-driven structures created at larger radii is ab- 
sorbed at the appropriate ILR. These changes to the dy- 
namical properties of the inner disk seem to be responsible 
for gradually increasing partial reflection of later incoming 
waves, and this partial feed-back in turn enhances the am- 
plitudes of subsequent disturbances. Growth is slow until 
the inner disk reflects waves strongly enough to create a 
coherent linear instability. 

Note that the slow, fluctuating rise in Smux in both runs 
50M and 500M (upper panel of Figure ^ is roughly ex- 
ponential, with similar exponents in both cases. This ap- 
pears to suggest that the behavior just described is, in 
fact, destabilizing in the sense that it implies that the rate 
of growth is independent of the amplitude. It would seem, 
therefore, that no system of particles, however large, could 
behave as a smooth disk. 

5. SUMMARY AND DISCUSSION 

The main result presented here is the demonstration 
that scattering at the inner Lindblad resonance of a spi- 
ral disturbance of even very low amplitude causes a lasting 
change to the properties of the disk that leads to increased 
amplitude of subsequent activity. Linear perturbation the- 
ory could never capture this behavior, because it neglects 
second order changes to the equilibrium model caused by 
a small amplitude disturbance. 

The random density fluctuations of a system of self- 
gravitating particles in a cool, shearing d isk are ampli- 
fied as they swing from le ading to trailing (lToomrelll990l : 
iToomre fc KalnaisI Il991l ) . The wave action created by 
these collective responses is carried inwards at the group 
velocity to the ILR where it is absorbed by particles. The 
resulting scattering of particles to more eccentric orbits 
de-populates the near-circular orbits over a narrow range 
of initial angular momentum. Here I have shown that, no 
matter how small the amplitude, the lasting changes to 
the particle distribution promote a higher level of density 
fluctuations that leads to indefinite growth. 

The growth of non-axisymmetric waves in the idealized 
simulations presented here is characterized by two phases. 
When very large numbers of particles are employed, the 
beha vior appears to be that of a systern of dr essed parti- 
cles (iToomre fc Kalnaislll99il : IWei nberslll99^. but one in 
which the amplitudes of the particle wakes rise slowly. As 
the peak relative overdensity passes ~ 2%, the changes to 
the background state are destabilizing and simple insta- 
bilities appear that exhibit run-away growth. 

The instability that appears is caused by earlier changes 
to the DF, and is a true mode that could have grown ex- 
ponentially from low amplitude. Direct tests reported in 
Section 14.31 demonstrate that it does not rely upon non- 
linear coupling to previous structures. I speculate on a 
possible mechanism for the un stable mode in Section [4.51 

Linear perturbation theory (|Toomrdll981l ) predicts no 
instability for a smooth DF, yet the evidence from Sec- 



tion 14.61 is that no finite number of particles would ever 
avoid indefinite growth. Whether or not this conclusion is 
correct, the avoidance of run-away growth for some non- 
infinite number of particles is of theoretical interest only, as 
galaxies are probably never as smooth as a disk of 5 x 10^ 
randomly distributed particles that are each only ^ 10 
times as massive as a typical star. 

Similar evolutionary changes to the original smooth disk 
are undoubtedly the cause of the non-axisymmetric fea- 
tu res that develope d for m > 2 in the simulations reported 
bylSellwood' ("2 01 If) to test the models proposed bv Bertin 
et al. (1989). The appearance of visible multi-arm waves 
that heated the disk in those models was, as here, more 
and more delayed as the number of particles was increased. 
Note that no large amplitude waves developed when forces 
were restricted to to = 2 because the disk was of such low 
mass. The swing amplification of shot noise is minimal 
when the parameter X > 3, which is the case for bisym- 
metric waves in this low-mass disk, and resonant scatter- 
ing by such weak disturbances will cause extremely slow 
growth, if any. However, X ex m~^ , making the disk more 
responsive to disturbances with m > 2, allowing run-away 
growth of multi-arm features and disk heating once these 
force terms were included. 

Toomre & Kalnajs (119911 ) used local simulations in a 
small shearing patch of a disk to study the behavior of 
swing-amplified shot noise. They were able to explain even 
the long-term behavior with just linear theory, although 
the limiting amplitude was a factor of ^ 2 higher than 
they estimated, which they attributed to additional corre- 
lations between the positions of particles that developed 
in the long-term. Most significantly, they did not observe 
the kind of run-away growth reported here. This difference 
may be due to a combination of factors: first, the particles 
in their principal box interacted with pair-wise forces and, 
had the box extended infinitely in the azimuthal direc- 
tion, each non-axisymmetric feature would be composed 
of a continuous spectrum of waves, implying that there 
would be no equivalent of Lindblad resonances and the 
absorption of wave action would be spread over a broad 
swath of particles. However, their periodic boundary con- 
ditions implied that the spectrum of waves was discrete, 
and a Lindblad resonance must have been present for each 
of the multiple waves that made up an individual shear- 
ing feature. Thus the absorption of wave action must have 
taken place at a number of discrete locations, and mild res- 
onant scattering probably occurred at each. But the most 
important distinguishing feature of their simulations that 
prevented development of destabilizing features at these 
resonances is that they applied a damping term to only 
the radial velocity component of their particles. Any res- 
onant scattering that depopulated nearly circular orbits 
must have been quickly erased, preventing their system 
from manifesting the behavior described here. 

I am grateful to the referee for sugges ting I note the ob - 
servations of Saturn's rings (reviewed in lCuzzi et aLll2010l ). 
which seem to support the picture of spiral chaos proposed 
by TbomrQ (1990 ). Of course, the snowball particles in 
this beautiful ring system are believed to collide dissipa- 
tively about once per orbit, which would again prevent 
scattering at possible Lindblad resonances from creating 
reflecting regions that might destabilize the disk. 
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While identifying the origin of indefinite growth is a 
good beginning to a detailed understanding the origin of 
spiral features in A^-body simulations, it is far from be- 
ing the whole story. All the simulations reported here 
supported several co-existing waves of differing frequen- 
cies at all times, except runs 50Mc and 50Md in which 
pre-existing waves were eliminated and single disturbance 
stood out for a brief period. The recurrence of large- 
amplitude waves remains to be described in a later paper 
of this series. 

In order to be able to reach the level of understanding 
described here, I have simplified the dynamics of a real- 
istic disk galaxy to the extreme. First, the particles in 
the simulations presented here interacted with each other 
only through the m = 2 sectoral harmonic of the distur- 
bance forces, with all other self- gravity terms being ig- 
nored; even the mean central attraction was replaced by 
that of a rigid mass distribution. Second, the particles 
were constrained to move in a plane only. Third, I delib- 
erately chose to stu dy an idealized model that was known 
(|Zana)ll976l : [Toomrelll98ll : lEvans fc Readlll998f ) to possess 
no small-amplitude instabilities. Fourth, I selected parti- 
cles so that their joint distribution of E and L^ matches 
that of the adopted DF as closely as possible. Fifth, the 
simulations were entirely isolated from external perturba- 
tions by passing satellites, halo substructure, etc., and the 
halo was represented by an idealized axisymmetric and 
time independent mass distribution. Sixth, all effects of 
gas dynamics, such as the creation of dense clouds and 
the formation of stars, were excluded. 

Yet even with all these abstractions, the behavior took 
some effort to understand. The demonstration of a sim- 
ple, clean instability was possible only in restarted simula- 
tions, constructed by careful reshuffling of the particles at 
the crucial time. Simulations that include all the above- 
mentioned complications continue to manifest recurrent 
spiral patterns. It is reasonable to hope that their ori- 
gin is related to the mechanism presented here, but this 
remains to be shown. 

The ultimate goal of this work is to demonstrate that 
some spirals in galaxies are self-excited i nstabilities with a 
simi lar origin to that i n the simulations. ISellwoodI (J2010b( ) 
and lHahn et al\ (J201l[ ) have made a step in this direction, 
by finding a feature in the distribution of solar neighbor- 
hood stars that resembles that of scattering at an ILR. Fur- 
ther analysis of data from other ground- and space-base d 
surveys, such as HERMES ("Bland-Hawthorn e t alll2010[ ). 
APOGEE (Maicwski et al. 2010) and especiall y Gaia C Per- 
ryman et al. |2001D . will afford better tests. 
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